Primary Outcome

Death from any cause or requirement of new intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support in the 28 days after randomisation.

This script was used to undertake the analysis of the primary outcome for ASCOT.

Preamble

Load packages

library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(plotly)
library(lubridate)
library(ggdist)
library(patchwork)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))
bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))
color_scheme_set("red")
options(digits = 4)

Prepare datasets

devtools::load_all()
# Raw data
all_dat <- read_all_no_daily()
all_dat_daily <- read_all_daily()

# FAS-ITT
fas_itt_dat <- all_dat %>% 
  filter_fas_itt() %>%
  transmute_model_cols_grp_aus_nz()
# FAS-ITT excluding missing values
fas_itt_nona_dat <- fas_itt_dat %>% 
  filter(!is.na(PO))

# ACS-ITT
acs_itt_dat <- all_dat %>% 
  filter_acs_itt() %>%
  transmute_model_cols_grp_aus_nz()
# ACS-ITT excluding missing values
acs_itt_nona_dat <- acs_itt_dat %>% 
  filter(!is.na(PO))
# ACS-ITT worst case
acs_itt_wc_dat <- acs_itt_dat %>%
  mutate(PO = if_else(is.na(PO) | PO == 1, 1, 0))
# ACS-ITT best case
acs_itt_bc_dat <- acs_itt_dat %>%
  mutate(PO = if_else(is.na(PO) | PO == 0, 0, 1))

# ACS-PP - exclude participants not per-protocol
pp <- read_per_protocol_list()
all_pp <- all_dat %>% left_join(pp, by = "StudyPatientID")
acs_pp_dat <- all_pp %>%
  filter_acs_itt() %>%
  filter(is.na(Reason)) %>%
  transmute_model_cols_grp_aus_nz()
acs_pp_nona_dat <- acs_pp_dat %>%
  filter(!is.na(PO))

# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
  filter_concurrent_intermediate()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>% 
  filter(!is.na(PO))
# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
  filter_concurrent_std_aspirin()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>% 
  filter(!is.na(PO))
# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
  filter_concurrent_therapeutic()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>% 
  filter(!is.na(PO))

Load models

logistic_mod <- cmdstan_model("stan/binary/logistic.stan")
logistic_site_mod <- cmdstan_model("stan/binary/logistic_site.stan")
logistic_epoch_mod <- cmdstan_model("stan/binary/logistic_epoch.stan")
logistic_site_epoch_llik <- cmdstan_model("stan/binary/logistic_site_epoch_loglik.stan")
logistic_site_epoch_llik_int <- cmdstan_model("stan/subgroup/logistic_site_epoch_loglik_shrinkage.stan")
primary_mod <- cmdstan_model("stan/binary/logistic_site_epoch.stan")

Helper functions

make_stan_data <- function(
    dat, 
    vars, 
    outcome, 
    beta_sd = c(2.5, rep(1, nXassign), 2.5, 10),
    ctr = contr.orthonorm) 
{
  X <- model.matrix(
     as.formula(paste("~", paste(vars, collapse = " + "))),
    data = dat,
    contrast = list(
      AAssignment = ctr,
      CAssignment = ctr)
  )
  nXassign <- sum(grepl("Assignment", colnames(X)))

  y <- dat[[outcome]]
  N <- dim(X)[1]
  K <- dim(X)[2]  
    epoch  <- dat$epoch
  M_epoch  <- max(dat$epoch)
    region <- dat[["ctry_num"]]
  M_region <- max(region)
    site <- dat[["site_num"]]
  M_site <- max(site)
  region_by_site <- region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)
  
  list(N = N, K = K, X = X, y = y,
       M_region = M_region, region = region,
       M_site = M_site, site = site,
       M_epoch = M_epoch, epoch = epoch,
       region_by_site = region_by_site,
       beta_sd = beta_sd)
}


# Always includes the assigned interventions
make_primary_model_data <- function(
    dat, 
    vars = c("inelgc3", "agegte60", "ctry"),
    beta_sd_var = c(10, 2.5, 1, 1),
    ctr = contr.orthonorm) {
  
  # Domain specific intercept as in interims
  XA <- make_domA_design(dat, ctr)
  # Domain specific intercept as in interims
  XC <- make_domC_design(dat, ctr)
  if (any(XC[, 1] != 1)) {
    add_intercept = TRUE
  } else {
    add_intercept = FALSE
  }
  
  if (!is.null(vars)) {
    Xother <- model.matrix(
     as.formula(paste(" ~ ", paste(vars, collapse = " + "))),
     data = dat
    ) [, -1]
  } else {
    Xother <- NULL
  }
  
  X <- cbind(XC, XA, Xother)
  if(add_intercept) {
    X <- cbind(randPlatform = 1, X)
  }
  nXtrt <- sum(grepl("rand", colnames(X))) - 1
  
  epoch  <- dat[["epoch"]]
  M_epoch  <- max(epoch)
  region <- dat[["ctry_num"]]
  M_region <- max(region)
  site <- dat[["site_num"]]
  M_site <- max(site)
  region_by_site <- region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)
  
  y <- dat[["PO"]]
  N <- dim(X)[1]
  K <- dim(X)[2]
  beta_sd <- c(2.5, rep(1, nXtrt), beta_sd_var)
  out <- list(
    N = N, K = K, X = X, y = y,
    M_region = M_region, region = region,
    M_site = M_site, site = site,
    M_epoch = M_epoch, epoch = epoch,
    region_by_site = region_by_site,
    beta_sd = beta_sd,
    ctrC = attr(XC, "contrasts")[[1]],
    ctrA = attr(XA, "contrasts")[[1]]
  )
  return(out)
}


fit_primary_model <- function(
    dat,
    model = primary_mod,
    vars = NULL,
    beta_sd_var = NULL,
    ctr = contr.orthonorm,
    seed = 32915,
    ...
) {
  mdat <- make_primary_model_data(
    dat,
    vars,
    beta_sd_var,
    ctr
  )
  snk <- capture.output(
    mfit <- model[["sample"]](
      data = mdat,
      seed = seed,
      refresh = 0,
      iter_warmup = 1000,
      iter_sampling = 2500,
      chains = 8,
      parallel_chains = min(8, parallel::detectCores()),
      ...
  ))
  mpars <- mfit$metadata()$model_params
  keep <- mpars[!grepl("(_raw|epsilon_)", mpars)]
  mdrws <- as_draws_rvars(mfit$draws(keep))
  names(mdrws$beta) <- colnames(mdat$X)
  
  if(any(grepl("site", keep))) {
    site_map <- dat %>% dplyr::count(site_num, site)  
    names(mdrws$gamma_site) <- site_map$site
  }
  if(any(grepl("epoch", keep))) {
    epoch_map <- dat %>% dplyr::count(epoch, epoch_lab)  
    names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
  }
  # Transformed samples
  Ca <- mdat$ctrA
  Cc <- mdat$ctrC
  # Get constrained parameters from contrast transformation
  mdrws$Acon <- Ca %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
  mdrws$Ccon <- as.vector(Cc %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))])
  names(mdrws$Ccon) <- rownames(Cc)
  # Get treatment effect in case of non treatment-coding
  mdrws$Atrt <- mdrws$Acon[-1] - mdrws$Acon[1]
  mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
  # Get odds ratios
  mdrws$OR <- exp(mdrws$Ctrt)
  mdrws$OR <- c(mdrws$OR, exp(mdrws$beta[!grepl("(Intercept|rand)", names(mdrws$beta))]))
  return(list(dat = mdat, fit = mfit, drws = mdrws))
}


make_po_table <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3", "C4"), 
      labels = intervention_labels()$CAssignment[-1])
    ) %>%
    group_by(CAssignment) %>%
    summarise(
      Randomised = n(),
      `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(PO)), 100 * mean(is.na(PO))),
      `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(PO)), 100 * mean(!is.na(PO))),
      `Met primary outcome` = sprintf("%i (%.1f)", sum(PO, na.rm = TRUE), 100 * mean(PO, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Randomised = n(),
          `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(PO)), 100 * mean(is.na(PO))),
          `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(PO)), 100 * mean(!is.na(PO))),
          `Met primary outcome` = sprintf("%i (%.1f)", sum(PO, na.rm = TRUE), 100 * mean(PO, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(CAssignment, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(
      font_size = 9, 
      latex_options = "HOLD_position")  
}


summarise_posterior <- function(dat, futval = 1.1) {
  dat %>%
    mutate(
      Median = 
        sprintf("%.2f", median(Posterior)),
      `95\\% CrI` = 
        sprintf("(%.2f, %.2f)", quantile(Posterior, 0.025), quantile(Posterior, 0.975)),
      `Mean (SD)` = 
        sprintf("%.2f (%.2f)", mean(Posterior), sd(Posterior)),
      `Pr(OR < 1)` = 
        sprintf("%.2f", Pr(Posterior < 1)),
      `Pr(OR > 1/1.1)` = 
        sprintf("%.2f", Pr(Posterior > 1/futval))
    )
}


make_subgroup_summary_table <- function(dat, format = "html") {
  L <- nrow(dat)
  row_st <- seq(1, L, by = L / 3)
  row_ed <- seq(L/3, L, by = L / 3)
  kable(
    dat[, -1] %>% select(-Posterior),
    format = format,
    booktabs = TRUE,
    align = "lrrrrr",
      escape = FALSE
  ) %>%
    kable_styling(
      latex_options = "HOLD_position",
      font_size = 9
    ) %>%
    group_rows("Intermediate", row_st[1], row_ed[1]) %>%
    pack_rows("Low with aspirin", row_st[2], row_ed[2]) %>%
    pack_rows("Therapeutic", row_st[3], row_ed[3])  
}

plot_or_densities <- function(rvs) {
  tibble(Contrast = fct_inorder(names(rvs)), RV = rvs) %>% 
  ggplot(., aes(y = Contrast, xdist = RV)) + 
  stat_halfeye(
    aes(fill = 
          stat(cut_cdf_qi(
            cdf, 
            .width = c(.5, .8, .95, 0.99), 
            labels = scales::percent_format()))),
    adjust = 1, n = 1001, .width = c(0.5, 0.8, 0.95)
    ) + 
  scale_fill_brewer(
    palette = "Reds",
    direction = -1, 
    na.translate = FALSE) + 
  labs(
    x = "Odds ratio contrast",
    fill = "Interval"
  ) +
  scale_x_log10(breaks = c(0.1, 0.25, 0.5, 1, 2, 4, 10)) + 
  geom_vline(xintercept = 1)
}


odds_ratio_summary_table <- function(OR, format = "html", fn = NULL) {
  out <- tibble(
    Parameter = names(OR),
    Median = median(OR),
    `95% CrI` = apply(
      quantile(OR, c(0.025, 0.975)), 2, 
      function(x) sprintf("(%.2f, %.2f)", x[1], x[2])),
    `Mean (SD)` = sprintf("%.2f (%.2f)", E(OR), sd(OR)),
    `Pr(OR < 1)` = Pr(OR < 1),
  ) %>%
  kable(
    format = format, 
    digits = 2, 
    align = "lrrrr",
    linesep = "",
    booktabs = TRUE) %>%
  kable_styling(
    font_size = 9,
    bootstrap_options = "striped", 
    latex_options = "HOLD_position")
  if (!is.null(fn) & format == "latex") {
    save_tex_table(out, fn)
  } else {
    return(out)
  }
}


decision_quantity_summary_table <- function(OR, format = "html", fut_val = 1.1) {
  tdat <- tibble(
    Intervention = names(OR),
    posterior = OR,
    Posterior = sprintf("%.2f (%.2f, %.2f)", 
                        median(posterior), 
                        quantile(posterior, 0.025), 
                        quantile(posterior, 0.975)),
    `Superior\nPr(OR = min(OR))` = sprintf("%.2f", E(posterior == rvar_min(posterior))),
    `Effective\nPr(OR < 1)` = sprintf("%.2f", Pr(posterior < 1)),
    `Futile\nPr(OR > 1/1.1)` = sprintf("%.2f", Pr(posterior > 1/fut_val)),
    `Equivalent\nPr(1/1.1 < OR < 1.1)` = sprintf("%.2f", Pr(posterior < fut_val & posterior > 1/fut_val))
  ) %>%
    select(-posterior)
  tdat[1, 2] <- "1.00"
  tdat[1, -(1:3)] <- "-"
  if(format == "latex") {
    colnames(tdat) <- linebreak(colnames(tdat), align = "c", linebreaker = "\n")  
  }
  out <- kable(
    tdat,
    format = format, 
    digits = 2, 
    align = "lrrrrr",
    linesep = "",
    escape = F,
    booktabs = TRUE) %>%
    kable_styling(
      font_size = 9,
      bootstrap_options = "striped",
      latex_options = "HOLD_position"
    )
  return(out)
}

Outline

This document presents the results of analyses of the primary outcome. The first section introduces the primary outcome definition and it’s derivation from the database fields. Descriptive analyses are then presented followed by the model-based analyses.

For the model-based analyses, results are presented for each population and sub-model.

Primary Outcome Definition

The primary outcome is a composite of death, need for new respiratory support, or vasopressor/inotropic support. From the protocol:

Death from any cause or requirement of new intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support in the 28 days after randomisation. This includes any participant who receives non-invasive mechanical ventilation (either CPAP or BIPAP, apart from the below considerations) any time after enrolment even if not transferred to ICU. It does NOT include the use of humidified high-flow nasal prong oxygen.

Participants on pre-existing home BiPAP or CPAP will not be considered to have met the primary outcome unless they have either:

• required invasive mechanical ventilation (i.e. intubation), or
• graduated from CPAP only whilst asleep to BiPAP at any time, or
• graduated from BiPAP only whilst asleep to BiPAP for >12 hours/day, or
• died by day 28

There may be cases where a patient has been assessed as requiring intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support, but the patient or family declined treatment and the patient was discharged home. If attempts to obtain 28-day data are unsuccessful or not possible, and the investigator had deemed at the time of discharge that the patient would be highly likely to die within 28 days from randomisation, these participants will be deemed to have met the primary outcome.

Derivation of the Outcome

Derivation of the outcome requires checking of the daily, discharge, and day 28 data extracts. On the daily data, there is a variable DD_PrimaryEndpointReachedToday however this was coded incorrectly in the original database and therefore fails to capture some participants. Additionally, given the composite nature of the outcome it is useful to check all components individually as well as the composite outcome. Therefore, this variable is not used in the derivation, but is included as a cross-check.

Below, each component is summarised in aggregate.

Day 28 mortality

Summarise the death component

mort_table <- all_dat %>%
  filter_acs_itt() %>%
  group_by(CAssignment) %>%
  summarise(
    `Mortality_Alive at day 28_` = sum(1 - D28_death, na.rm = TRUE),
    `Mortality_Death within 28 days_Total` = sum(D28_death, na.rm = TRUE),
    `Mortality_Death within 28 days_Prior to discharge` = sum(DIS_death & DIS_day <= 28, na.rm = TRUE),
    `Mortality_Death within 28 days_Post-discharge` = sum(D28_death == 1 & DIS_death == 0, na.rm = TRUE),
    `Mortality_Unknown_` = sum(is.na(D28_death))
  ) %>%
  gather(key, Frequency, -CAssignment, factor_key = T) %>%
  spread(CAssignment, Frequency) %>%
  separate(key, into = c("Measure", "Outcome", "Breakdown"), sep = "_") %>%
  mutate(
    Overall = C1 + C2 + C3 + C4,
    across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x[c(1,2,5)])))
  )
kable(
  mort_table %>% select(-1),
  caption = "Composite primary endpoint breakdown - mortality.",
  align = "llrrrrr") %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Mortality", 1, nrow(mort_table))

Outcome	Breakdown	C1	C2	C3	C4	Overall
Mortality
Alive at day 28		577 (95)	588 (96)	271 (96)	44 (88)	1480 (95)
Death within 28 days	Total	19 (3)	15 (2)	10 (4)	6 (12)	50 (3)
Death within 28 days	Prior to discharge	15 (2)	11 (2)	10 (4)	4 (8)	40 (3)
Death within 28 days	Post-discharge	4 (1)	4 (1)	0 (0)	2 (4)	10 (1)
Unknown		14 (2)	10 (2)	2 (1)	0 (0)	26 (2)

Composite primary endpoint breakdown - mortality.

Vasopressor/inotropic requirement

Summarise the vasopressor/inotropes component

vaso_table <- all_dat %>%
  filter_acs_itt() %>%
  group_by(CAssignment) %>%
  summarise(
    `Vasopressor/inotropes_Not required_` = sum(1 - ANY_vasop, na.rm = TRUE),
    `Vasopressor/inotropes_Use within 28 days_Total` = sum(ANY_vasop, na.rm = TRUE),
    `Vasopressor/inotropes_Use within 28 days_Prior to discharge` = sum(DD_vasop, na.rm = TRUE),
    `Vasopressor/inotropes_Use within 28 days_Post-discharge` = sum(D28_vasop, na.rm = TRUE),
    `Vasopressor/inotropes_Unknown_` = sum(is.na(ANY_vasop))
  ) %>%
  gather(key, Frequency, -CAssignment, factor_key = T) %>%
  spread(CAssignment, Frequency) %>%
  separate(key, into = c("Measure", "Outcome", "Breakdown"), sep = "_") %>%
  mutate(
    Overall = C1 + C2 + C3 + C4,
    across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x[c(1,2,5)])))
  )
kable(
  vaso_table %>% select(-1),
  caption = "Composite primary endpoint breakdown - vasopressor/inotropes.",
  align = "llrrrr") %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Vasopressor/inotropes", 1, nrow(vaso_table))

Outcome	Breakdown	C1	C2	C3	C4	Overall
Vasopressor/inotropes
Not required		591 (97)	595 (97)	273 (96)	48 (96)	1507 (97)
Use within 28 days	Total	6 (1)	6 (1)	6 (2)	2 (4)	20 (1)
Use within 28 days	Prior to discharge	5 (1)	6 (1)	6 (2)	2 (4)	19 (1)
Use within 28 days	Post-discharge	2 (0)	0 (0)	0 (0)	1 (2)	3 (0)
Unknown		13 (2)	12 (2)	4 (1)	0 (0)	29 (2)

Composite primary endpoint breakdown - vasopressor/inotropes.

New intensive respiratory support

Summarise the ventilation component

vent_table <- all_dat %>%
  filter_acs_itt()  %>%
  group_by(CAssignment) %>%
  summarise(
    `Ventilation_Not required_` = sum(1 - ANY_vent, na.rm = TRUE),
    `Ventilation_Use within 28 days_Total` = sum(ANY_vent, na.rm = TRUE),
    `Ventilation_Use within 28 days_Prior to discharge` = sum(ANY_DD_vent, na.rm = TRUE),
    `Ventilation_Use within 28 days_Post-discharge` = sum(ANY_vent & ANY_DD_vent == 0, na.rm = TRUE),
    `Ventilation_Unknown_` = sum(is.na(ANY_vent))
  ) %>%
  gather(key, Frequency, -CAssignment, factor_key = T) %>%
  spread(CAssignment, Frequency) %>%
  separate(key, into = c("Measure", "Outcome", "Breakdown"), sep = "_") %>%
  mutate(
    Overall = C1 + C2 + C3 + C4,
    across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x[c(1, 2, 5)])))
  )
kable(
  vent_table %>% select(-1),
  caption = "Composite primary endpoint breakdown - respiratory support.",
  align = "llrr") %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Ventilation", 1, nrow(vent_table))

Outcome	Breakdown	C1	C2	C3	C4	Overall
Ventilation
Not required		562 (92)	579 (94)	263 (93)	43 (86)	1447 (93)
Use within 28 days	Total	35 (6)	24 (4)	18 (6)	7 (14)	84 (5)
Use within 28 days	Prior to discharge	29 (5)	19 (3)	13 (5)	5 (10)	66 (4)
Use within 28 days	Post-discharge	6 (1)	5 (1)	5 (2)	2 (4)	18 (1)
Unknown		13 (2)	10 (2)	2 (1)	0 (0)	25 (2)

Composite primary endpoint breakdown - respiratory support.

DAMA and Likely to Die

There were 3 participants who were DAMA and identified as likely to die. However, the day 28 status was observed for all 3 participants, therefore, no-one met the primary endpoint due to DAMA and unknown day 28 status.

Summarise the DAMA component

all_dat %>%
  filter_acs_itt()  %>%
  select(DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
  dplyr::count(DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
  spread(D28_death, n) %>%
  kable(
    col.names = c("DAMA", "Likely to die by day 28", "No", "Yes", "(Missing)")
  ) %>%
  kable_styling() %>%
  add_header_above(c(" " = 2, "Death by day 28" = 3))

		Death by day 28
DAMA	Likely to die by day 28	No	Yes	(Missing)
0	0	1447	46	24
1	0	32	2	2
1	1	1	2	NA

Composite

Summarise the composite

comp_table <- all_dat %>%
  filter_acs_itt()  %>%
  group_by(CAssignment) %>%
  summarise(
    `Primary Outcome_No_`  = sum(1 - PO, na.rm = TRUE),
    `Primary Outcome_Yes_` = sum(PO, na.rm = TRUE),
    `Primary Outcome_Unknown_Total` = sum(is.na(PO)),
    `Primary Outcome_Unknown_Day 28 status` = sum(is.na(PO) & is.na(D28_death)),
    `Primary Outcome_Unknown_Vasopressor/inotropes` = 
      sum(is.na(PO) & !is.na(D28_death) & is.na(ANY_vasop))
  ) %>%
  gather(key, Frequency, -CAssignment, factor_key = T) %>%
  spread(CAssignment, Frequency) %>%
  separate(key, into = c("Measure", "Outcome", "Breakdown"), sep = "_") %>%
  mutate(
    Overall = C1 + C2 + C3 + C4,
    across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x[1:3])))
  )
kable(
  comp_table %>% select(-1),
  align = "llrr") %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Primary outcome", 1, nrow(comp_table))

Outcome	Breakdown	C1	C2	C3	C4	Overall
Primary outcome
No		561 (92)	576 (94)	259 (92)	43 (86)	1439 (92)
Yes		35 (6)	25 (4)	20 (7)	7 (14)	87 (6)
Unknown	Total	14 (2)	12 (2)	4 (1)	0 (0)	30 (2)
Unknown	Day 28 status	13 (2)	10 (2)	2 (1)	0 (0)	25 (2)
Unknown	Vasopressor/inotropes	1 (0)	2 (0)	2 (1)	0 (0)	5 (0)

Table 1: Composite primary endpoint breakdown - overall.

Summarise the all components

tab <- bind_rows(comp_table, mort_table, vaso_table, vent_table)
tab_head <- all_dat %>% 
  filter_acs_itt() %>% 
  dplyr::count(CAssignment = factor(CAssignment, labels = intervention_labels()$CAssignment[-1])) %>%
  mutate(label = paste0(CAssignment, "<br>(n = ", n, ")")) %>%
  pull(label) %>%
  c(paste0("Overall<br><br>(n = ", nrow(all_dat %>% filter_acs_itt()), ")"))
colnames(tab)[4:8] <- tab_head
kable(
  tab %>% select(-1),
  align = "llrrrrr",
  escape = F) %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Primary outcome", 1, nrow(comp_table)) %>%
  group_rows("Mortality", nrow(comp_table) + 1, nrow(comp_table) + nrow(mort_table)) %>%
  group_rows("Vasopressor/inotropes", 11, 15) %>%
  group_rows("Ventilation", 16, 20)

Outcome	Breakdown	Low dose (n = 610)	Intermediate dose (n = 613)	Low dose with aspirin (n = 283)	Therapeutic dose (n = 50)	Overall (n = 1556)
Primary outcome
No		561 (92)	576 (94)	259 (92)	43 (86)	1439 (92)
Yes		35 (6)	25 (4)	20 (7)	7 (14)	87 (6)
Unknown	Total	14 (2)	12 (2)	4 (1)	0 (0)	30 (2)
Unknown	Day 28 status	13 (2)	10 (2)	2 (1)	0 (0)	25 (2)
Unknown	Vasopressor/inotropes	1 (0)	2 (0)	2 (1)	0 (0)	5 (0)
Mortality
Alive at day 28		577 (95)	588 (96)	271 (96)	44 (88)	1480 (95)
Death within 28 days	Total	19 (3)	15 (2)	10 (4)	6 (12)	50 (3)
Death within 28 days	Prior to discharge	15 (2)	11 (2)	10 (4)	4 (8)	40 (3)
Death within 28 days	Post-discharge	4 (1)	4 (1)	0 (0)	2 (4)	10 (1)
Unknown		14 (2)	10 (2)	2 (1)	0 (0)	26 (2)
Vasopressor/inotropes
Not required		591 (97)	595 (97)	273 (96)	48 (96)	1507 (97)
Use within 28 days	Total	6 (1)	6 (1)	6 (2)	2 (4)	20 (1)
Use within 28 days	Prior to discharge	5 (1)	6 (1)	6 (2)	2 (4)	19 (1)
Use within 28 days	Post-discharge	2 (0)	0 (0)	0 (0)	1 (2)	3 (0)
Unknown		13 (2)	12 (2)	4 (1)	0 (0)	29 (2)
Ventilation
Not required		562 (92)	579 (94)	263 (93)	43 (86)	1447 (93)
Use within 28 days	Total	35 (6)	24 (4)	18 (6)	7 (14)	84 (5)
Use within 28 days	Prior to discharge	29 (5)	19 (3)	13 (5)	5 (10)	66 (4)
Use within 28 days	Post-discharge	6 (1)	5 (1)	5 (2)	2 (4)	18 (1)
Unknown		13 (2)	10 (2)	2 (1)	0 (0)	25 (2)

Table 2: Composite primary endpoint breakdown.

Save table to outputs

colnames(tab)[4:8] <- linebreak(tab_head, align = "c", linebreaker = "<br>")
out <- kable(
  tab %>% select(-1),
  format = "latex",
  align = "llrrrrr",
  booktabs = TRUE,
  escape = F) %>%
  kable_styling(font_size = 9, latex_options = "HOLD_position") %>%
  collapse_rows(1, valign = "top", latex_hline = 'none') %>%
  group_rows("Primary outcome", 1, nrow(comp_table)) %>%
  group_rows("Mortality", nrow(comp_table) + 1, nrow(comp_table) + nrow(mort_table)) %>%
  group_rows("Vasopressor/inotropes", 11, 15) %>%
  group_rows("Ventilation", 16, 20)
save_tex_table(out, fn = "outcomes/primary/composite-breakdown")

Summarise the component combinations

levs <- c("None", "Death + Vas./Ino. + Ventilation", "Death + Vas./Ino.", "Death + Ventilation", "Vas./Ino. + Ventilation", 
          "Death", "Vas./Ino.", "Ventilation", "Unknown")
tab <- all_dat %>%
  filter_acs_itt()  %>%
  group_by(CAssignment) %>%
  summarise(
    `Composite outcomes` = case_when(
      PO == 0 ~ 0,
      is.na(PO) ~ 99,
      D28_death == 1 & ANY_vasop == 1 & ANY_vent == 1 ~ 1,
      D28_death == 1 & ANY_vasop == 1 & ANY_vent == 0 ~ 2,
      D28_death == 1 & ANY_vasop == 0 & ANY_vent == 1 ~ 3,
      D28_death == 0 & ANY_vasop == 1 & ANY_vent == 1 ~ 4,
      D28_death == 1 ~ 5,
      ANY_vasop == 1 ~ 6,
      ANY_vent == 1 ~ 7,
    )) %>%
  dplyr::count(CAssignment, `Composite outcomes`) %>%
  mutate(`Composite outcomes` = factor(`Composite outcomes`, levels = c(0:7, 99), labels = levs)) %>%
  spread(CAssignment, n, fill = 0) %>%
  mutate(Overall = C1 + C2 + C3 + C4) %>%
  mutate(across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x))))
colnames(tab)[-1] <- tab_head
colnames(tab)[1] <- "Composite outcomes, n (\\%)"
kable(
  tab ,
  align = "lrrrrr",
  escape = F) %>%
  kable_styling(font_size = 11)

Composite outcomes, n (\%)	Low dose (n = 610)	Intermediate dose (n = 613)	Low dose with aspirin (n = 283)	Therapeutic dose (n = 50)	Overall (n = 1556)
None	561 (92)	576 (94)	259 (92)	43 (86)	1439 (92)
Death + Vas./Ino. + Ventilation	6 (1)	3 (0)	3 (1)	2 (4)	14 (1)
Death + Ventilation	13 (2)	11 (2)	7 (2)	4 (8)	35 (2)
Vas./Ino. + Ventilation	0 (0)	2 (0)	1 (0)	0 (0)	3 (0)
Death	0 (0)	1 (0)	0 (0)	0 (0)	1 (0)
Vas./Ino.	0 (0)	1 (0)	2 (1)	0 (0)	3 (0)
Ventilation	16 (3)	7 (1)	7 (2)	1 (2)	31 (2)
Unknown	14 (2)	12 (2)	4 (1)	0 (0)	30 (2)

Table 3: Observed outcome combinations.

Code

colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
out <- kable(
  tab,
  format = "latex",
  align = "lrrrrr",
  booktabs = TRUE,
  linesep = "",
  escape = F) %>%
  kable_styling(
    font_size = 9, 
    latex_options = "HOLD_position")
save_tex_table(out, fn = "outcomes/primary/composite-combinations")

Descriptive Analyses

By Interim Analysis

Code

intdat <- all_dat %>%
  filter_acs_itt() %>%
  mutate(
    interim = as.character(case_when(
      RandDate < get_interim_dates()$meet_date[1] ~ 1,
      RandDate < get_interim_dates()$meet_date[2] ~ 2,
      RandDate < get_interim_dates()$meet_date[3] ~ 3,
      RandDate < get_interim_dates()$meet_date[4] ~ 4,
      TRUE ~ 4
    )),
    CAssignment = factor(
      CAssignment, labels = c("Low", "Intermediate", "Low\nwith aspirin", "Therapuetic"))
  )
intcounts <- intdat %>% 
  dplyr::count(interim) %>% 
  mutate(label = paste0("Interim ", interim, " (n = ", n, ")"))
trtcounts <- intdat %>% dplyr::count(CAssignment)
ovrtab <- intdat %>%
  group_by(interim = "Overall", CAssignment) %>%
  summarise(
    Randomised = 
      sprintf("%i", n()),
    Known = 
      sprintf("%i (%.0f)", sum(!is.na(PO)), 100 * mean(!is.na(PO))),
    `Met primary outcome` = 
      sprintf("%i (%.0f)", sum(PO, na.rm = T), 100 * mean(PO, na.rm = T)),
    `Death` = 
      sprintf("%i (%.0f)", sum(D28_death, na.rm = T), 100 * mean(D28_death, na.rm = T)),
    `Vasopressor/inotropic support` = 
      sprintf("%i (%.0f)", sum(ANY_vasop, na.rm = T), 100 * mean(ANY_vasop, na.rm = T)),
    `Intensive respiratory support` = 
      sprintf("%i (%.0f)", sum(ANY_vent, na.rm = T), 100 * mean(ANY_vent, na.rm = T))
  ) %>%
  pivot_longer(`Randomised`:`Intensive respiratory support`, names_to = "Outcome") %>%
  pivot_wider(names_from = CAssignment, values_from = value, values_fill = "-")
inttab <- intdat %>%
  group_by(interim, CAssignment) %>%
  summarise(
    Randomised = 
      sprintf("%i", n()),
    Known = 
      sprintf("%i (%.0f)", sum(!is.na(PO)), 100 * mean(!is.na(PO))),
    `Met primary outcome` = 
      sprintf("%i (%.0f)", sum(PO, na.rm = T), 100 * mean(PO, na.rm = T)),
    `Death` = 
      sprintf("%i (%.0f)", sum(D28_death, na.rm = T), 100 * mean(D28_death, na.rm = T)),
    `Vasopressor/inotropic support` = 
      sprintf("%i (%.0f)", sum(ANY_vasop, na.rm = T), 100 * mean(ANY_vasop, na.rm = T)),
    `Intensive respiratory support` = 
      sprintf("%i (%.0f)", sum(ANY_vent, na.rm = T), 100 * mean(ANY_vent, na.rm = T))
  ) %>%
  pivot_longer(`Randomised`:`Intensive respiratory support`, names_to = "Outcome") %>%
  pivot_wider(names_from = CAssignment, values_from = value, values_fill = "-")
tab <- bind_rows(ovrtab, inttab)

savetab <- kable(
  tab[, -1],
  format = "latex",
  align = "lrrrr",
  booktabs = TRUE
) %>%
  kable_styling(
    font_size = 9,
    latex_options = "HOLD_position"
  ) %>%
  group_rows(paste0("Overall (n = ", sum(intcounts$n), ")"), 1, 6) %>%
  group_rows(intcounts$label[1], 7, 12) %>%
  group_rows(intcounts$label[2], 13, 18) %>%
  group_rows(intcounts$label[3], 19, 24) %>%
  group_rows(intcounts$label[4], 25, 30)
save_tex_table(savetab, fn = "outcomes/primary/primary-composite-interims")

kable(
  tab[, -1],
  format = "html",
  align = "lrrrr",
) %>%
  kable_styling(
    font_size = 9,
    latex_options = "HOLD_position"
  ) %>%
  group_rows(paste0("Overall (n = ", sum(intcounts$n), ")"), 1, 6) %>%
  group_rows(intcounts$label[1], 7, 12) %>%
  group_rows(intcounts$label[2], 13, 18) %>%
  group_rows(intcounts$label[3], 19, 24) %>%
  group_rows(intcounts$label[4], 25, 30)

Outcome	Low	Intermediate	Low with aspirin	Therapuetic
Overall (n = 1556)
Randomised	610	613	283	50
Known	596 (98)	601 (98)	279 (99)	50 (100)
Met primary outcome	35 (6)	25 (4)	20 (7)	7 (14)
Death	19 (3)	15 (2)	10 (4)	6 (12)
Vasopressor/inotropic support	6 (1)	6 (1)	6 (2)	2 (4)
Intensive respiratory support	35 (6)	24 (4)	18 (6)	7 (14)
Interim 1 (n = 685)
Randomised	234	226	225	-
Known	227 (97)	222 (98)	222 (99)	-
Met primary outcome	22 (10)	15 (7)	19 (9)	-
Death	13 (6)	11 (5)	9 (4)	-
Vasopressor/inotropic support	2 (1)	3 (1)	6 (3)	-
Intensive respiratory support	22 (10)	15 (7)	17 (8)	-
Interim 2 (n = 246)
Randomised	90	98	58	-
Known	90 (100)	97 (99)	57 (98)	-
Met primary outcome	2 (2)	3 (3)	1 (2)	-
Death	1 (1)	3 (3)	1 (2)	-
Vasopressor/inotropic support	1 (1)	0 (0)	0 (0)	-
Intensive respiratory support	2 (2)	3 (3)	1 (2)	-
Interim 3 (n = 384)
Randomised	177	172	-	35
Known	173 (98)	167 (97)	-	35 (100)
Met primary outcome	8 (5)	4 (2)	-	6 (17)
Death	4 (2)	0 (0)	-	5 (14)
Vasopressor/inotropic support	3 (2)	2 (1)	-	2 (6)
Intensive respiratory support	8 (5)	4 (2)	-	6 (17)
Interim 4 (n = 241)
Randomised	109	117	-	15
Known	106 (97)	115 (98)	-	15 (100)
Met primary outcome	3 (3)	3 (3)	-	1 (7)
Death	1 (1)	1 (1)	-	1 (7)
Vasopressor/inotropic support	0 (0)	1 (1)	-	0 (0)
Intensive respiratory support	3 (3)	2 (2)	-	1 (7)

Age

Relationship age to outcome

agedat <- acs_itt_dat %>%
  dplyr::count(PO, AgeAtEntry) %>% 
  spread(PO, n, fill = 0) %>% 
  mutate(
    n = `0` + `1` + `<NA>`,
    p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p1 <- ggplot(agedat, aes(AgeAtEntry, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry")
p2 <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_vline(xintercept = 60, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\nwith outcome", x = "Age at entry")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-age.pdf"), p, height = 2, width = 6)
p

Figure 1: Relationship (logistic regression linear in age) between age at entry and the primary outcome.

Relationship age 265 60 to outcome

acs_itt_dat %>%
  dplyr::count(`Age 60+` = agegte60, PO) %>%
  group_by(`Age 60+`) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  kable(digits = 2) %>%
  kable_styling(bootstrap_options = "striped")

Age 60+	0	1	<NA>	p_1	p_miss
0	1034	47	21	0.04	0.02
1	405	40	9	0.09	0.02

Table 4: Relationship between age >= 60 and the primary outcome.

Country

Relationship country to outcome

tdat <- all_dat %>%
  filter_acs_itt() %>%
  dplyr::count(Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")), PO) %>%
  group_by(Country) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(tdat, aes(Country, p_1)) +
  geom_point() +
    labs(
      y = "Proportion meeting\nprimary outcome", 
      x = "Country of enrolment")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-country.pdf"), p, height = 2, width = 6)
p

Figure 2: Primary outcome by country.

Site

Relationship site to outcome

tdat <- all_dat %>%
  filter_acs_itt() %>%
  dplyr::count(
    Country_lab = Country,
    Site_lab = fct_infreq(Location),
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")),
    Site = PT_LocationName,
    PO) %>%
  group_by(Country, Site) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()

p1 <- ggplot(tdat, aes(Site_lab, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(tdat, aes(Site_lab, p_1)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_point() +
    labs(
      y = "Proportion meeting\nprimary outcome", 
      x = "Site of enrolment") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p <- p1 / p2
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-country-site.pdf"), p, height = 4, width = 6.25)
p

Figure 3: Primary outcome by site within country.

Calendar Time

Relationship calendar date to outcome

caldat <- acs_itt_dat %>%
  dplyr::count(PO, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(PO, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`),
         n = `1` + `0` + `<NA>`)
p1 <- ggplot(caldat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_point() +
    labs(
      y = "Proportion meeting\nprimary outcome", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12) +
  ylim(0, 0.25)
p2 <- ggplot(caldat, aes(mth, n)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p <- p2 | p1
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-calendar-time.pdf"), p, height = 2, width = 6)
p

Figure 4: Relationship between calendar date and the primary outcome.

Days since symptom onset

Relationship days since sympton onset to outcome

dsfsdat <- acs_itt_dat %>%
  dplyr::count(PO, dsfs) %>% 
  spread(PO, n, fill = 0) %>% 
  mutate(n = `1` + `0` + `<NA>`,
         p = `1` / (`1` + `0`))
dsfsmod <- glm(
  cbind(`1`, `0`) ~ dsfs, 
  data = dsfsdat, 
  family = binomial())
dsfsdat <- dsfsdat %>%
  mutate(
    ypred = predict(dsfsmod, newdata = dsfsdat, type = "response")
  )
p1 <- ggplot(dsfsdat, aes(dsfs, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 7.5, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Days since first symptoms")
p2 <- ggplot(dsfsdat, aes(dsfs, p)) +
    geom_point() +
    geom_vline(xintercept = 7.5, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\nwith outcome", 
         x = "Days since first symptoms")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-dsfs.pdf"), p, height = 2, width = 6)
p

Figure 5: Relationship (logistic regression linear) between days since symptom onset at entry and the primary outcome.

Anti-coagulation Domain

Relationship anti-coagulation to outcome

make_po_table(
    all_dat %>%
      filter_acs_itt()
)

n (\%)	Low dose	Intermediate dose	Low dose with aspirin	Therapeutic dose	Overall
Randomised	610	613	283	50	1556
Outcome missing	14 (2.3)	12 (2.0)	4 (1.4)	0 (0.0)	30 (1.9)
Outcome observed	596 (97.7)	601 (98.0)	279 (98.6)	50 (100.0)	1526 (98.1)
Met primary outcome	35 (5.9)	25 (4.2)	20 (7.2)	7 (14.0)	87 (5.7)

Table 5: Primary outcome by anti-coagulation intervention.

Relationship anti-coagulation to outcome

save_tex_table(make_po_table(
      all_dat %>%
      filter_acs_itt(),
      "latex"), 
      "outcomes/primary/anticoagulation-summary")

Analyses

ACS-ITT

In these analyses, the population is restricted to participants who were randomised to the anti-coagulation domain.

Treatment only

Data and sampling treatment only model

res <- fit_primary_model(acs_itt_nona_dat, logistic_mod, seed = 17250)
names(res$drws$OR) <- intervention_labels2()$CAssignment[-(1:2)]

Code

plot_or_densities(res$drws$OR)

Figure 6: Posterior densities for odds ratio contrasts.

Code

odds_ratio_summary_table(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.71	(0.43, 1.20)	0.74 (0.20)	0.90
Low-dose with aspirin	1.23	(0.70, 2.14)	1.28 (0.37)	0.23
Therapeutic-dose	2.40	(0.95, 5.34)	2.60 (1.16)	0.03

Table 6: Posterior summaries for odds ratio contrasts.

Diagnostics

Code

res$fit$summary("beta")

# A tibble: 6 × 10
  variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
  <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -2.62   -2.61   0.143 0.143 -2.86   -2.39    1.00   15499.   14290.
2 beta[2]   0.485   0.494  0.334 0.330 -0.0778  1.02    1.00   17188.   14447.
3 beta[3]  -0.210  -0.210  0.210 0.211 -0.554   0.136   1.00   17721.   15370.
4 beta[4]  -0.695  -0.694  0.246 0.245 -1.10   -0.289   1.00   17656.   15552.
5 beta[5]  -0.0510 -0.0333 0.399 0.396 -0.735   0.580   1.00   20107.   14608.
6 beta[6]  -0.785  -0.774  0.504 0.501 -1.63    0.0235  1.00   21997.   13319.

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0843 1.0378 1.0546 1.1228 1.0945 1.0608 0.9765 1.0343

Trace plots

Code

mcmc_trace(res$drws["beta"])

Treatment plus covariates (no site or epoch)

Under this model a component is added for country (region), \(\mathbb{I}[\text{age}\geq 60]\), and intervention specific ineligibility (C3 only)

New terms:

\[ \begin{aligned} \rho_1 &= 0 \\ \rho_r &\overset{\text{iid}}{\sim} \text{Normal}(0, 1), r=2,...,R,\\ \beta_{\text{age}\geq 60} &\sim \text{Normal}(0, 2.5^2) \\ \beta_{\text{inelgC3}} &\sim \text{Normal}(0, 10) \end{aligned} \]

Code

res <- fit_primary_model(
  acs_itt_nona_dat, 
  logistic_mod, 
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 71236)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Code

plot_or_densities(res$drws$OR)

Figure 7: Posterior densities for odds ratio contrasts.

Code

odds_ratio_summary_table(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.71	(0.42, 1.21)	0.74 (0.20)	0.90
Low-dose with aspirin	1.32	(0.73, 2.33)	1.37 (0.41)	0.18
Therapeutic-dose	1.85	(0.68, 4.60)	2.05 (1.03)	0.11
Ineligible aspirin	2.39	(0.74, 6.29)	2.67 (1.44)	0.07
Age ≥ 60	1.96	(1.25, 3.05)	2.01 (0.46)	0.00
Australia/New Zealand	1.00	(0.29, 3.06)	1.17 (0.74)	0.50
Nepal	1.43	(0.65, 2.92)	1.53 (0.59)	0.18

Table 7: Posterior summaries for odds ratio contrasts.

Diagnostics

Code

res$fit$summary("beta")

# A tibble: 10 × 10
   variable    mean   median    sd   mad      q5      q95  rhat ess_bulk
   <chr>      <dbl>    <dbl> <dbl> <dbl>   <dbl>    <dbl> <dbl>    <dbl>
 1 beta[1]  -3.00   -3.00    0.198 0.197 -3.33   -2.68     1.00   14028.
 2 beta[2]   0.253   0.255   0.378 0.378 -0.379   0.866    1.00   15820.
 3 beta[3]  -0.155  -0.158   0.216 0.216 -0.506   0.198    1.00   19065.
 4 beta[4]  -0.626  -0.628   0.261 0.262 -1.05   -0.193    1.00   17111.
 5 beta[5]   0.0869  0.0802  0.646 0.645 -0.960   1.15     1.00   16336.
 6 beta[6]  -0.795  -0.787   0.497 0.495 -1.62    0.00501  1.00   21338.
 7 beta[7]   0.843   0.873   0.542 0.532 -0.0981  1.68     1.00   22546.
 8 beta[8]   0.674   0.675   0.226 0.224  0.298   1.05     1.00   18977.
 9 beta[9]  -0.0176 -0.00389 0.606 0.600 -1.03    0.953    1.00   16486.
10 beta[10]  0.352   0.358   0.382 0.381 -0.290   0.965    1.00   18103.
# … with 1 more variable: ess_tail <dbl>

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0215 1.0098 1.0926 1.0081 1.0257 0.9988 1.0248 1.0574

Trace plots

Code

mcmc_trace(res$drws["beta"])

Treatment plus covariates plus site (no epoch)

Under this model a component is added for site nested within country.

New terms

\[ \begin{aligned} \xi_{rs} &\overset{\text{iid}}{\sim}\text{Normal}(0,\sigma^2_{r}), s=1,...,S_r, r=1,...,R \\ \sigma_r &\overset{\text{iid}}{\sim} \text{Half-}t(3, 1). \end{aligned} \]

Data and sampling no site or epoch model

res <- fit_primary_model(
  acs_itt_nona_dat, 
  logistic_site_mod, 
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 52765,
  adapt_delta = 0.95)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Code

plot_or_densities(res$drws$OR)

Figure 8: Posterior densities for odds ratio contrasts.

Code

odds_ratio_summary_table(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.72	(0.41, 1.24)	0.75 (0.21)	0.88
Low-dose with aspirin	0.94	(0.51, 1.74)	0.99 (0.31)	0.57
Therapeutic-dose	2.07	(0.73, 5.66)	2.37 (1.31)	0.08
Ineligible aspirin	2.72	(0.75, 8.52)	3.20 (2.07)	0.06
Age ≥ 60	1.80	(1.11, 2.93)	1.86 (0.46)	0.01
Australia/New Zealand	1.01	(0.22, 4.21)	1.32 (1.10)	0.49
Nepal	1.87	(0.46, 6.63)	2.30 (1.70)	0.17

Table 8: Posterior summaries for odds ratio contrasts.

Code

plot_site_terms(
  exp(res$drws$gamma_site), 
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)

Figure 9: Site effect odds ratios.

Diagnostics

Code

res$fit$summary(c("beta", "gamma_site", "tau_site"))

# A tibble: 39 × 10
   variable     mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>       <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -3.56    -3.55   0.532 0.519 -4.45   -2.71    1.00    5496.    8830.
 2 beta[2]   0.576    0.577  0.404 0.399 -0.101   1.23    1.00   17006.   14724.
 3 beta[3]  -0.0975  -0.0956 0.228 0.228 -0.474   0.280   1.00   18674.   15698.
 4 beta[4]  -0.526   -0.526  0.272 0.272 -0.974  -0.0808  1.00   17773.   14426.
 5 beta[5]   0.195    0.196  0.700 0.702 -0.961   1.35    1.00   16795.   14547.
 6 beta[6]  -0.966   -0.955  0.547 0.542 -1.88   -0.0897  1.00   23767.   14794.
 7 beta[7]   0.978    1.00   0.619 0.613 -0.0588  1.95    1.00   22493.   14369.
 8 beta[8]   0.588    0.586  0.246 0.243  0.184   0.992   1.00   25540.   15554.
 9 beta[9]   0.00242  0.0117 0.750 0.745 -1.26    1.20    1.00   12505.   14538.
10 beta[10]  0.614    0.627  0.667 0.645 -0.510   1.68    1.00    9105.   12175.
# … with 29 more rows

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 2

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8537 0.8557 0.8084 0.8207 0.7391 0.8056 0.8420 0.7953

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["tau_site"])

Treatment plus covariates plus epoch (no site)

Replace site terms with epoch terms

\[ \begin{aligned} \tau_1 &= 0 \\ \tau_t &= \tau_{t-1} + \sigma_\tau\epsilon_t \\ \epsilon_t &\overset{\text{iid}}{\sim} \text{Normal}(0,1) \\ \sigma_\tau &\sim \text{Half-}t(3, 1) \end{aligned} \]

Data and sampling no site model

res <- fit_primary_model(
  acs_itt_nona_dat, 
  logistic_epoch_mod, 
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 76405,
  adapt_delta = 0.98)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Code

plot_or_densities(res$drws$OR)

Figure 10: Posterior densities for odds ratio contrasts.

Code

odds_ratio_summary_table(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.73	(0.43, 1.24)	0.76 (0.21)	0.87
Low-dose with aspirin	0.92	(0.50, 1.68)	0.97 (0.31)	0.60
Therapeutic-dose	2.30	(0.83, 6.08)	2.59 (1.39)	0.05
Ineligible aspirin	1.70	(0.52, 4.59)	1.92 (1.07)	0.17
Age ≥ 60	1.95	(1.22, 3.09)	2.00 (0.48)	0.00
Australia/New Zealand	1.41	(0.40, 4.58)	1.70 (1.11)	0.29
Nepal	2.91	(1.15, 7.19)	3.23 (1.58)	0.01

Table 9: Posterior summaries for odds ratio contrasts.

Code

plot_epoch_terms(
  exp(res$drws$gamma_epoch)
)

Figure 11: Epoch effect odds ratios.

Diagnostics

Code

res$fit$summary(c("beta", "gamma_epoch", "tau_epoch"))

# A tibble: 23 × 10
   variable   mean median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>     <dbl>  <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.22  -4.19  0.528 0.511 -5.13   -3.39    1.00   15486.   12840.
 2 beta[2]   0.659  0.662 0.406 0.402 -0.0127  1.32    1.00   18941.   15955.
 3 beta[3]  -0.125 -0.122 0.219 0.219 -0.486   0.235   1.00   20213.   15097.
 4 beta[4]  -0.535 -0.534 0.266 0.267 -0.976  -0.102   1.00   18884.   15031.
 5 beta[5]   0.348  0.340 0.654 0.651 -0.709   1.43    1.00   16689.   14620.
 6 beta[6]  -0.833 -0.819 0.512 0.507 -1.70   -0.0150  1.00   25915.   14359.
 7 beta[7]   0.504  0.530 0.553 0.547 -0.431   1.37    1.00   24106.   13640.
 8 beta[8]   0.663  0.667 0.237 0.236  0.273   1.05    1.00   25646.   14940.
 9 beta[9]   0.339  0.346 0.622 0.625 -0.707   1.34    1.00   17486.   14074.
10 beta[10]  1.06   1.07  0.465 0.464  0.291   1.82    1.00   17849.   14963.
# … with 13 more rows

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8896 0.7985 0.8545 0.8386 0.8910 0.8580 0.8231 0.8728

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_epoch"])

Code

mcmc_trace(res$drws["tau_epoch"])

Primary Model

Data and sampling primary model

res <- fit_primary_model(
  acs_itt_nona_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 32915,
  adapt_delta = 0.99
)
res$fit$save_object(file.path("outputs", "models", "primary", "acs-itt-primary.rds"))
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Code

save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-summary-table")
odds_ratio_summary_table(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.74	(0.43, 1.27)	0.77 (0.22)	0.86
Low-dose with aspirin	0.88	(0.47, 1.64)	0.93 (0.30)	0.65
Therapeutic-dose	2.22	(0.77, 6.20)	2.55 (1.45)	0.07
Ineligible aspirin	2.50	(0.69, 7.93)	2.96 (1.94)	0.08
Age ≥ 60	1.88	(1.16, 3.04)	1.94 (0.49)	0.01
Australia/New Zealand	1.05	(0.23, 4.39)	1.37 (1.13)	0.47
Nepal	2.17	(0.51, 7.67)	2.66 (1.97)	0.13

Posterior summaries for model parameters.

Code

save_tex_table(
  decision_quantity_summary_table(c("Low-dose" = rvar(1), res$drws$OR[1:3]), "latex"),
  "outcomes/primary/primary-model-acs-itt-decision-table")
decision_quantity_summary_table(
  c("Low-dose" = rvar(1), res$drws$OR[1:3])
)

Intervention	Posterior	Superior Pr(OR = min(OR))	Effective Pr(OR < 1)	Futile Pr(OR > 1/1.1)	Equivalent Pr(1/1.1 < OR < 1.1)
Low-dose	1.00	0.08	-	-	-
Intermediate-dose	0.74 (0.43, 1.27)	0.63	0.86	0.23	0.15
Low-dose with aspirin	0.88 (0.47, 1.64)	0.27	0.65	0.47	0.22
Therapeutic-dose	2.22 (0.77, 6.20)	0.02	0.07	0.95	0.04

Table 10: Decision quantity summaries

Code

p <- plot_or_densities(res$drws$OR[1:3]) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "primary-model-acs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Figure 12: Posterior densities for odds ratio contrasts.

Code

p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "primary", "primary-model-acs-itt-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Diagnostics

Code

res$fit$summary(c("beta", "gamma_epoch", "tau_epoch", "gamma_site", "tau_site"))

# A tibble: 52 × 10
   variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.07   -4.06   0.671 0.654 -5.21   -3.01    1.00    8389.   12597.
 2 beta[2]   0.666   0.668  0.420 0.416 -0.0262  1.36    1.00   20492.   15814.
 3 beta[3]  -0.105  -0.103  0.228 0.227 -0.483   0.270   1.00   22714.   15494.
 4 beta[4]  -0.502  -0.499  0.277 0.277 -0.961  -0.0501  1.00   20392.   14539.
 5 beta[5]   0.224   0.218  0.698 0.697 -0.910   1.37    1.00   19342.   15216.
 6 beta[6]  -0.947  -0.944  0.542 0.533 -1.84   -0.0635  1.00   27617.   15229.
 7 beta[7]   0.899   0.915  0.624 0.618 -0.162   1.89    1.00   27309.   14051.
 8 beta[8]   0.634   0.633  0.248 0.249  0.222   1.04    1.00   27924.   14818.
 9 beta[9]   0.0392  0.0515 0.755 0.752 -1.24    1.26    1.00   14945.   14640.
10 beta[10]  0.749   0.773  0.687 0.683 -0.428   1.83    1.00   11783.   13803.
# … with 42 more rows

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8537 0.7930 0.7605 0.7244 0.7954 0.7820 0.7752 0.7430

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_epoch"])

Code

mcmc_trace(res$drws["tau_epoch"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["tau_site"])

Posterior Predictive Checks

Code

y_ppc <- res$drws$y_ppc
ppc_dat <- bind_cols(acs_itt_nona_dat, tibble(y_ppc = y_ppc))

grp_ppc <- function(grp) {
  ppc_dat %>%
    group_by(grp = !!grp) %>%
    summarise(
      mean_y = mean(PO),
      rvar_mean_y_ppc = rvar_mean(y_ppc)
    )
}
plot_grp_ppc <- function(dat, lab = "") {
  dat %>%
    ggplot(aes(y = grp, xdist = rvar_mean_y_ppc)) +
    stat_interval(size = 2) +
    geom_point(aes(x = mean_y, y = 1:nrow(dat)), colour = "red", shape = 23) +
    labs(
      x = "Posterior predictive\nproportion", 
      y = lab)  
}

ppc_C <- grp_ppc(sym("CAssignment"))
ppc_ctry <- grp_ppc(sym("country"))
ppc_epoch <- grp_ppc(sym("epoch"))
ppc_site <- ppc_dat %>%
  group_by(Country = country, grp = site_raw) %>%
  summarise(
    mean_y = mean(PO),
    rvar_mean_y_ppc = rvar_mean(y_ppc)
  )

p1 <- plot_grp_ppc(ppc_C, "Anticoagulation\nintervention") + labs(x = "")
p2 <- plot_grp_ppc(ppc_ctry, "Country") + labs(x = "")
p3 <- plot_grp_ppc(ppc_epoch, "Epoch") + labs(x = "")
p4 <- plot_grp_ppc(ppc_site %>% filter(Country == "IN"), "Sites\nIndia")  + labs(x = "")
p5 <- plot_grp_ppc(ppc_site %>% filter(Country == "AU"), "Sites\nAustralia") + labs(x = "")
p6 <- plot_grp_ppc(ppc_site %>% filter(Country == "NP"), "Sites\nNepal")
p7 <- plot_grp_ppc(ppc_site %>% filter(Country == "NZ"), "Sites\nNZ")
p <- ((p3 | p1 / p2) / 
        ( (p4 | p5) / (p6 | p7) + 
            plot_layout(heights = c(3, 1)) ) ) +
  plot_layout(heights = c(1, 1.5), guides = "collect")
pth <- file.path("outputs", "figures", "outcomes", "primary",
                 "primary-model-acs-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.5)
p

Sensitivity Analyses

Sens: Treatment Contrasts

Code

res <- fit_primary_model(
  acs_itt_nona_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  ctr = contr.treatment,
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Code

odds_ratio_summary_table(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.71	(0.42, 1.21)	0.74 (0.20)	0.90
Low-dose with aspirin	0.85	(0.46, 1.55)	0.89 (0.28)	0.70
Therapeutic-dose	2.04	(0.72, 5.39)	2.29 (1.22)	0.08
Ineligible aspirin	2.44	(0.66, 7.74)	2.89 (1.92)	0.08
Age ≥ 60	1.90	(1.17, 3.07)	1.96 (0.49)	0.01
Australia/New Zealand	1.18	(0.25, 4.84)	1.52 (1.27)	0.41
Nepal	2.10	(0.50, 7.44)	2.58 (1.93)	0.14

Posterior summaries for model parameters under treatment contrast priors.

Code

save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-trt-ctr-summary-table")

Code

p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "primary",
                 "primary-model-acs-itt-trt-ctr-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Sens: Weakly Informative Priors

As a prior sensitivity assessment, the primary model is re-fit under weakly informative priors (Normal(0, 10)) for all fixed-effect parameters.

Code

res <- fit_primary_model(
  acs_itt_nona_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 10, 10, 10),
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Code

odds_ratio_summary_table(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.74	(0.43, 1.29)	0.77 (0.22)	0.86
Low-dose with aspirin	0.89	(0.47, 1.65)	0.93 (0.30)	0.64
Therapeutic-dose	2.10	(0.73, 5.83)	2.41 (1.36)	0.08
Ineligible aspirin	2.51	(0.67, 8.19)	2.99 (2.00)	0.08
Age ≥ 60	1.89	(1.16, 3.10)	1.95 (0.50)	0.01
Australia/New Zealand	1.37	(0.09, 12.45)	2.59 (4.28)	0.40
Nepal	4.05	(0.54, 32.02)	7.29 (17.05)	0.07

Posterior summaries for model parameters under treatment contrast priors.

Code

p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
p

Sens: FAS-ITT

In these analyses, the population is includes all participants who were randomised to the platform, even if not randomised into the anticoagulation domain.

Code

res <- fit_primary_model(
  fas_itt_nona_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Figure 13: Posterior densities for odds ratio contrasts.

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.74	(0.43, 1.27)	0.77 (0.22)	0.86
Low-dose with aspirin	0.89	(0.48, 1.62)	0.93 (0.29)	0.64
Therapeutic-dose	2.17	(0.77, 6.00)	2.48 (1.37)	0.07
Ineligible aspirin	2.16	(0.60, 6.75)	2.55 (1.64)	0.11
Age ≥ 60	1.89	(1.18, 3.05)	1.95 (0.48)	0.01
Australia/New Zealand	1.11	(0.26, 4.56)	1.43 (1.16)	0.44
Nepal	2.19	(0.50, 7.89)	2.70 (2.00)	0.13

Table 11: Posterior summaries for odds ratio contrasts.

Code

save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-fas-itt-summary-table")

Code

p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "primary", "primary-model-fas-itt-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Sens: Worst-Case

Under this scenario, all participants with unknown primary outcome status were were to assumed to have met the primary outcome.

Code

make_po_table(acs_itt_wc_dat)

n (\%)	Low dose	Intermediate dose	Low dose with aspirin	Therapeutic dose	Overall
Randomised	610	613	283	50	1556
Outcome missing	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)
Outcome observed	610 (100.0)	613 (100.0)	283 (100.0)	50 (100.0)	1556 (100.0)
Met primary outcome	49 (8.0)	37 (6.0)	24 (8.5)	7 (14.0)	117 (7.5)

Code

save_tex_table(
  make_po_table(acs_itt_wc_dat, "latex"), 
  "outcomes/primary/acs-itt-worst-case-anticoagulation-summary")

Code

res <- fit_primary_model(
  acs_itt_wc_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Posterior densities for odds ratio contrasts.

Code

save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-worst-case-summary-table")
odds_ratio_summary_table(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.74	(0.46, 1.19)	0.77 (0.19)	0.89
Low-dose with aspirin	0.81	(0.46, 1.40)	0.84 (0.24)	0.77
Therapeutic-dose	1.37	(0.50, 3.51)	1.53 (0.78)	0.26
Ineligible aspirin	1.22	(0.35, 3.67)	1.42 (0.89)	0.37
Age ≥ 60	1.67	(1.07, 2.59)	1.71 (0.39)	0.01
Australia/New Zealand	2.37	(0.54, 9.18)	2.99 (2.34)	0.12
Nepal	1.89	(0.43, 6.83)	2.31 (1.72)	0.18

Posterior summaries for model parameters.

Code

p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal")))
pth <- file.path("outputs", "figures", "outcomes", "primary", 
                 "primary-model-acs-itt-worst-case-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Sens: Best-Case

Under this scenario, all participants with unknown primary outcome status were were to assumed to have not met the primary outcome.

Code

save_tex_table(
  make_po_table(acs_itt_bc_dat, "latex"), 
  "outcomes/primary/acs-itt-best-case-anticoagulation-summary")
make_po_table(acs_itt_bc_dat)

n (\%)	Low dose	Intermediate dose	Low dose with aspirin	Therapeutic dose	Overall
Randomised	610	613	283	50	1556
Outcome missing	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)	0 (0.0)
Outcome observed	610 (100.0)	613 (100.0)	283 (100.0)	50 (100.0)	1556 (100.0)
Met primary outcome	35 (5.7)	25 (4.1)	20 (7.1)	7 (14.0)	87 (5.6)

Code

res <- fit_primary_model(
  acs_itt_bc_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Posterior densities for odds ratio contrasts.

Code

save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-best-case-summary-table")
odds_ratio_summary_table(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.75	(0.43, 1.30)	0.78 (0.22)	0.85
Low-dose with aspirin	0.91	(0.49, 1.67)	0.95 (0.31)	0.62
Therapeutic-dose	2.32	(0.80, 6.45)	2.66 (1.49)	0.06
Ineligible aspirin	2.83	(0.80, 9.03)	3.37 (2.19)	0.05
Age ≥ 60	1.86	(1.13, 3.04)	1.92 (0.49)	0.01
Australia/New Zealand	1.00	(0.22, 4.07)	1.29 (1.06)	0.50
Nepal	2.15	(0.49, 7.71)	2.64 (1.95)	0.14

Posterior summaries for model parameters.

Code

p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal")))
pth <- file.path("outputs", "figures", "outcomes", "primary", 
                 "primary-model-acs-itt-best-case-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Sens: Concurrent Randomisation

Intermediate-dose

Code

save_tex_table(
  make_po_table(
    acs_itt_concurc2_dat,
    "latex"), 
  "outcomes/primary/acs-itt-c2-summary")
make_po_table(acs_itt_concurc2_dat)

n (\%)	Low dose	Intermediate dose	Overall
Randomised	610	613	1223
Outcome missing	14 (2.3)	12 (2.0)	26 (2.1)
Outcome observed	596 (97.7)	601 (98.0)	1197 (97.9)
Met primary outcome	35 (5.9)	25 (4.2)	60 (5.0)

Table 12: Primary outcome summary for participants randomised to either low-dose or intermediate dose arms (ACS-ITT-intermediate).

Code

ctr <- contr.orthonorm
X <- model.matrix(
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc2_nona_dat,
  contrasts.arg = list(CAssignment = ctr))
y <- acs_itt_concurc2_nona_dat$PO
sdat <- list(
  N = nrow(X),
  K = ncol(X),
  X = X,
  y = y,
  beta_sd = c(2.5, 1, 2.5, 1, 1)
)
snk <- capture.output(
  sfit <- logistic_mod$sample(
    data = sdat,
    seed = 3274,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500
  )
)
sdrws <- as_draws_rvars(sfit$draws("beta"))
names(sdrws$beta) <- colnames(X)
sdrws$beta_trt <- as.vector(
  transform_coding(cbind(1, ctr(2)), cbind(1, contr.treatment(2))) %**% sdrws$beta[1:2]
)
sdrws$OR <- c(
  "Intermediate-dose" = exp(sdrws$beta_trt[2]),
  "Age \u2265 60" = exp(sdrws$beta[3]),
  "Australia/New Zealand" = exp(sdrws$beta[4]),
  "Nepal" = exp(sdrws$beta[5])
)

Code

save_tex_table(
  odds_ratio_summary_table(sdrws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-c2-summary-table")
odds_ratio_summary_table(sdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.69	(0.41, 1.16)	0.71 (0.19)	0.92
Age ≥ 60	2.19	(1.29, 3.66)	2.26 (0.61)	0.00
Australia/New Zealand	1.39	(0.64, 2.77)	1.47 (0.55)	0.20
Nepal	1.06	(0.40, 2.38)	1.15 (0.51)	0.45

Table 13: Parameter posterior summaries.

Diagnostics

Code

res$fit$summary(c("beta"))

# A tibble: 10 × 10
   variable    mean   median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>    <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.08   -4.07    0.667 0.652 -5.20   -3.01    1.00    8163.   11086.
 2 beta[2]   0.678   0.683   0.420 0.417 -0.0206  1.37    1.00   22134.   15754.
 3 beta[3]  -0.130  -0.131   0.229 0.227 -0.508   0.245   1.00   22749.   15590.
 4 beta[4]  -0.517  -0.517   0.275 0.276 -0.969  -0.0673  1.00   21721.   16386.
 5 beta[5]   0.258   0.261   0.698 0.698 -0.898   1.41    1.00   21338.   14916.
 6 beta[6]  -0.979  -0.968   0.542 0.542 -1.89   -0.112   1.00   30272.   14314.
 7 beta[7]   1.03    1.04    0.619 0.616 -0.0154  2.02    1.00   29722.   14530.
 8 beta[8]   0.620   0.620   0.252 0.252  0.204   1.03    1.00   32053.   14952.
 9 beta[9]  -0.0154 -0.00183 0.743 0.748 -1.26    1.18    1.00   14921.   15151.
10 beta[10]  0.742   0.767   0.692 0.669 -0.432   1.83    1.00   11785.   14311.

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7842 0.7528 0.7704 0.8127 0.7986 0.7521 0.7614 0.7696

Trace plots

Code

mcmc_trace(res$drws["beta"])

Low-dose with aspirin

Code

save_tex_table(
  make_po_table(
    acs_itt_concurc3_dat,
    "latex"), 
  "outcomes/primary/acs-itt-c3-summary")
make_po_table(acs_itt_concurc3_dat)

n (\%)	Low dose	Intermediate dose	Low dose with aspirin	Overall
Randomised	299	298	283	880
Outcome missing	7 (2.3)	5 (1.7)	4 (1.4)	16 (1.8)
Outcome observed	292 (97.7)	293 (98.3)	279 (98.6)	864 (98.2)
Met primary outcome	21 (7.2)	16 (5.5)	20 (7.2)	57 (6.6)

Table 14: Primary outcome summary for participants randomised whilst the low-dose plus aspirin arm was available (ACS-ITT-aspirin).

Code

X <- model.matrix(
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc3_nona_dat,
  contrasts.arg = list(CAssignment = contr.orthonorm))
y <- acs_itt_concurc3_nona_dat$PO
sdat <- list(
  N = nrow(X),
  K = ncol(X),
  X = X,
  y = y,
  beta_sd = c(2.5, 1, 1, 2.5, 1)
)
snk <- capture.output(
  sfit <- logistic_mod$sample(
    data = sdat,
    seed = 3274,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500
  )
)
sdrws <- as_draws_rvars(sfit$draws("beta"))
names(sdrws$beta) <- colnames(X)

sdrws$beta_trt <- as.vector(
  transform_coding(cbind(1, contr.orthonorm(3)), cbind(1, contr.treatment(3))) %**% sdrws$beta[1:3]
)

sdrws$OR <- c(
  "Intermediate-dose" = exp(sdrws$beta_trt[2]),
  "Low-dose with aspirin" = exp(sdrws$beta_trt[3]),
  "Age \u2265 60" = exp(sdrws$beta[4]),
  "Australia/New Zealand" = exp(sdrws$beta[5])
)

Code

save_tex_table(
  odds_ratio_summary_table(sdrws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-c3-summary-table")
odds_ratio_summary_table(sdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate-dose	0.75	(0.39, 1.45)	0.79 (0.28)	0.80
Low-dose with aspirin	0.93	(0.50, 1.73)	0.98 (0.32)	0.59
Age ≥ 60	1.91	(1.09, 3.28)	1.98 (0.57)	0.01
Australia/New Zealand	0.31	(0.07, 1.09)	0.38 (0.27)	0.97

Table 15: Parameter posterior summaries.

Therapeutic-dose

Code

save_tex_table(
  make_po_table(
    acs_itt_concurc4_dat,
    "latex"), 
  "outcomes/primary/acs-itt-c4-summary")
make_po_table(acs_itt_concurc4_dat)

n (\%)	Low dose	Intermediate dose	Therapeutic dose	Overall
Randomised	79	65	50	194
Outcome missing	4 (5.1)	3 (4.6)	0 (0.0)	7 (3.6)
Outcome observed	75 (94.9)	62 (95.4)	50 (100.0)	187 (96.4)
Met primary outcome	6 (8.0)	3 (4.8)	7 (14.0)	16 (8.6)

Table 16: Primary outcome summary for participants randomised under protocol 5.0 (ACS-ITT-therapeutic).

Code

X <- model.matrix(
  ~ CAssignment + agegte60 + ctry,
  data = acs_itt_concurc4_nona_dat,
  contrasts.arg = list(CAssignment = contr.orthonorm))
y <- acs_itt_concurc4_nona_dat$PO
sdat <- list(
  N = nrow(X),
  K = ncol(X),
  X = X,
  y = y,
  beta_sd = c(2.5, 1, 1, 2.5, 1, 1)
)
snk <- capture.output(
  sfit <- logistic_mod$sample(
    data = sdat,
    seed = 3274,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500
  )
)
sdrws <- as_draws_rvars(sfit$draws("beta"))
names(sdrws$beta) <- colnames(X)

sdrws$beta_trt <- as.vector(
  transform_coding(cbind(1, contr.orthonorm(3)), cbind(1, contr.treatment(3))) %**% sdrws$beta[1:3]
)

sdrws$OR <- c(
  "Intermediate dose" = exp(sdrws$beta_trt[2]),
  "Therapeutic dose" = exp(sdrws$beta_trt[3]),
  "Age \u2265 60" = exp(sdrws$beta[4]),
  "India" = exp(sdrws$beta[5]),
  "Australia/New Zealand" = exp(sdrws$beta[6])
)

Code

save_tex_table(
  odds_ratio_summary_table(sdrws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-c4-summary-table")
odds_ratio_summary_table(sdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate dose	0.65	(0.18, 2.25)	0.80 (0.57)	0.75
Therapeutic dose	1.63	(0.54, 5.07)	1.93 (1.20)	0.19
Age ≥ 60	1.88	(0.66, 5.45)	2.18 (1.27)	0.12
India	0.39	(0.08, 1.62)	0.51 (0.41)	0.90
Australia/New Zealand	1.05	(0.37, 2.84)	1.19 (0.65)	0.46

Table 17: Parameter posterior summaries.

Per-Protocol Analysis

In the per-protocol analysis, participants who did not satisfy the intervention protocol (as determined by blinded review) are excluded from the analysis set.

Code

all_pp %>%
  filter(ENR_rec == 1) %>%
  dplyr::count(Reason)

# A tibble: 4 × 2
  Reason                              n
  <chr>                           <int>
1 Patient Withdrawal from Domain     26
2 Protocol Deviation                 28
3 Protocol Deviation & Withdrawal     1
4 <NA>                             1544

Code

acs_itt_dat %>%
  left_join(pp, by = "StudyPatientID") %>%
  dplyr::count(CAssignment, Reason) %>%
  spread(CAssignment, n, fill = 0)

# A tibble: 4 × 5
  Reason                             C1    C2    C3    C4
  <chr>                           <dbl> <dbl> <dbl> <dbl>
1 Patient Withdrawal from Domain      2     4     1     0
2 Protocol Deviation                  9     6     8     3
3 Protocol Deviation & Withdrawal     0     0     0     1
4 <NA>                              599   603   274    46

Code

acs_itt_dat %>%
  left_join(pp, by = "StudyPatientID") %>%
  dplyr::count(country, Reason) %>%
  spread(country, n, fill = 0)

# A tibble: 4 × 5
  Reason                             IN    AU    NP    NZ
  <chr>                           <dbl> <dbl> <dbl> <dbl>
1 Patient Withdrawal from Domain      4     3     0     0
2 Protocol Deviation                 19     7     0     0
3 Protocol Deviation & Withdrawal     0     1     0     0
4 <NA>                             1265   115   118    24

Code

save_tex_table(
  make_po_table(
    acs_pp_dat,
    "latex"), 
  "outcomes/primary/acs-pp-outcome-summary")
make_po_table(acs_pp_dat)

n (\%)	Low dose	Intermediate dose	Low dose with aspirin	Therapeutic dose	Overall
Randomised	599	603	274	46	1522
Outcome missing	13 (2.2)	11 (1.8)	3 (1.1)	0 (0.0)	27 (1.8)
Outcome observed	586 (97.8)	592 (98.2)	271 (98.9)	46 (100.0)	1495 (98.2)
Met primary outcome	35 (6.0)	25 (4.2)	19 (7.0)	7 (15.2)	86 (5.8)

Table 18: Primary outcome summary for ACS-PP set.

Data and sampling primary model

mdat <- make_primary_model_data(
  acs_pp_nona_dat,
  c("inelgc3", "agegte60", "ctry"),
  c(10, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- primary_mod$sample(
    data = mdat,
    seed = 81241,
    refresh = 0,
    iter_warmup = 1000,
    iter_sampling = 2500,
    chains = 8,
    parallel_chains = 8,
    adapt_delta = 0.99
))
mdrws <- as_draws_rvars(mfit$draws(c(
  "beta", 
  "gamma_site", "tau_site", 
  "gamma_epoch", "tau_epoch", "y_ppc")))
names(mdrws$beta) <- colnames(mdat$X)

# Transformed samples
epoch_map <- acs_pp_nona_dat %>% dplyr::count(epoch, epoch_lab)
site_map <- acs_pp_nona_dat %>% dplyr::count(site_num, site)
site_facet <- factor(mdat$region_by_site, labels = c("India", "Australia\nNew Zealand", "Nepal"))

Ca <- mdat$ctrA
Cc <- mdat$ctrC

mdrws$Acon <- Ca %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
mdrws$Ccon <- Cc %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
mdrws$Atrt <- mdrws$Acon[-1] - mdrws$Acon[1]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]

mdrws$OR <- setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic"))
mdrws$OR <- c(
  mdrws$OR, 
  "Ineligible aspirin" = exp(mdrws$beta[7]), 
  "Age \u2265 60" = exp(mdrws$beta[8]),
   setNames(exp(mdrws$beta[9:10]), c("Australia/New Zealand", "Nepal")))
mdrws$ORepoch <- setNames(exp(mdrws$gamma_epoch), epoch_map$epoch_lab)
mdrws$ORsites <- setNames(exp(mdrws$gamma_site), site_map$site)

Code

save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/primary/primary-model-acs-pp-summary-table")
odds_ratio_summary_table(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Intermediate	0.74	(0.42, 1.28)	0.77 (0.22)	0.86
Low with aspirin	0.85	(0.44, 1.56)	0.89 (0.29)	0.70
Therapeutic	2.52	(0.86, 7.13)	2.91 (1.67)	0.05
Ineligible aspirin	2.46	(0.66, 8.05)	2.94 (1.98)	0.09
Age ≥ 60	1.87	(1.14, 3.05)	1.93 (0.49)	0.01
Australia/New Zealand	1.11	(0.24, 4.49)	1.43 (1.18)	0.44
Nepal	1.99	(0.49, 7.13)	2.47 (1.87)	0.16

Table 19: Posterior summaries for model parameters.

Code

p <- plot_epoch_site_terms(mdrws$ORepoch, mdrws$ORsites, site_facet)
pth <- file.path("outputs", "figures", "outcomes", "primary", 
                 "primary-model-acs-pp-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p